October 16 2024
The lab counts for up to 4 points towards the final grade of the course.
This lab assumes that some of the procedures learned in previous labs for downloading and manipulating images in R are already mastered and therefore they are not explained in detail here. Refer to previous labs if there is any problem with those procedures.
The study area corresponds to the vicinity of the city of Pucallpa. It is crossed by the Ucayali river which along with the Marañon river, form eventually the Amazon river. Pucallpa is connected to Lima, the capital of Peru, through the Federico Basadre road. This road triggered rapid deforestation for the establishment of different cultivations. Most of the agricultural activity has consisted of slash and burn cropping and most recently of the establishment of large and small oil palm plantations.
Place the image in your working directory. Then, open a new script in R Studio and change the route to your working directory in the script using the function setwd().
Copy and paste each chunk of code in your new R script and run it trying to understand the purpose, logic and syntaxis of each line. Make sure the code runs with no errors before moving to the next one.
Answer the questions in the answer sheet and submit it along with the required pdf file to canvas.
Have fun!
library(terra)
library(RStoolbox)
wd="/Users/tug61163/Library/CloudStorage/OneDrive-TempleUniversity/Courses/IntroRemoteSensing/2024Fall/Class7_Indices/LabMaterials"
setwd(wd)
getwd() # tells you what working directory you are in, just to be sure it's correct
tars=list.files(wd, pattern="tar")
tars
untar(tars)
TIFS=list.files(wd, pattern=".TIF")
TIFS
# Make sure the numbers corresponding to the names of bands 1 to 7 in the tifs object
L8 = rast(TIFS[c(3:9)]) # rast() function brings the raster bands into the R environment as rasters. Specifically, as a spatRaster object.
L8
plotRGB(L8, r=5, g=4, b=3, stretch="lin") # this may have to be run in the console directly. It must plot in the RStudio plot window, usually bottom right
e=draw(x = "extent") # click top left, bottom right
L8rsz=crop(L8, e)
plotRGB(L8rsz, r=5, g=4, b=3, stretch="lin")
# Remove objects to release memory
rm(L8)
# SPECTRAL INDICES
# Let's calculate the NDVI (normalized difference vegetation index) manually
# NDVI = (NIR - Red) / (NIR + Red)
NDVI = ((L8rsz[[5]] - L8rsz[[4]]) / (L8rsz[[5]] + L8rsz[[4]]))
plot(NDVI)
# You can use the function spectralIndices() to calculate a variety of indices. Check ?spectralIndices. We will illustrate with NDVI and two additional indices:
L8indices=spectralIndices(L8rsz, blue=2, green=3, red=4, nir=5, swir2=6, swir3=7, indices=c("NDVI", "NDWI", "NBRI"))
plot(L8indices[[1]])
plot(L8indices[[2]])
plot(L8indices[[3]])
plotRGB(L8indices, r=3, g=1, b=2, stretch="lin")
# Tasseled Cap requires to drop the first band in Landsat 8. See ?tasseledCap
L8tascap=tasseledCap(L8rsz[[-1]], sat="Landsat8OLI")
plotRGB(L8tascap, r=1, g=2, b=3, stretch="lin")
L8PCA=rasterPCA(L8rsz)
L8PCA
summary(L8PCA$model)
plotRGB(L8PCA$map, r=1, g=2, b=3, stretch="lin")
writeRaster(L8rsz, filename="L8rsz.tif")
writeRaster(L8indices, filename="L8indices.tif")
writeRaster(L8tascap, filename="L8tascap.tif")
writeRaster(L8PCA$map, filename="L8PCA.tif")
In the QGIS browser panel, right click on “XYZ Tiles”. Then select “New connections”.
This will open a new XYZ connection window. Enter “google satellite” as the Name and the following string as URL: https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}
This will add a new layer with high resolution imagery from google.
You can alternatively install the quickmap services plugin following the procedure that will be shown in class.
Zoom in to different areas representing these land covers: Urban Old-growth forest Lake River Grass Burnt areas Sand banks Oil palm
Then explore and interpret the colors representing each land cover in
the RGB maps produced with the spectral indices and the tasseled cap
transformation. The image below can be helpful to evaluate the
contribution of different bands to pixel color in RGB images:
Download from Earth Explorer a recent Landsat image covering the study area that you have considered for your final project. The image should have low cloud cover and good quality. If your study area is subjected to seasons, select an image from nearly the peak of the growing season (July-August). Crop the image to the size of your study area.
Produce a tasseled cap transformation and import the image to QGIS.
Use the imported image along with high-resolution google imagery to identify and list below the main land cover types that yo can discriminate in your study area.
Zoom in to a representative area for each land cover type and produce two screenshots per land cover:
What land cover results in the redest colors? ___________________________ Interpretation: ___________________________________________________________________
What land cover results in the bluest colors? ____________________________ Interpretation: __________________________________________________________________
What land cover results in the greenest color? _________________________________ Interpretation: __________________________________________________
are there any land covers appearing as yellow, cyan, magenta or white? which are those?_____________________________________________ Interpretation: _________________________________________________
BONUS QUESTION (0.5 e points).Upload upload an RGB map with bands 1, 2, and 3 resulting from a PCA transformation. Describe to what extent areas representing each one of the land cover classes that you defined for your study area tend to display consistently similar color characteristics .
To upload in canvas:
The name of the land cover types defined for the study area as described in point 3 (0.5 pt)
A single pdf file with the two screenshots for areas tipifying the land cover classes that you identified for your study area described above. Make sure to label the screenshots with the land cover type that they represent and to use an adequate zoom to allow appropriate visualization (2 pt).
Answer to the questions in point 5 (1.5 pt) in a word document.
OPTIONAL: answer to bonus question (0.5 pt)